Modeling of Interstellar Scintillation Arcs from Pulsar B1133+16 
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ABSTRACT 

The parabolic arc phenomenon visible in the Fourier analysis of the scintillation spectra of 
pulsars provides a new method of investigating the small scale structure in the ionized interstellar 
medium (ISM). We report archival observations of the pulsar B1133+16 showing both forward and 
reverse parabolic arcs sampled over 14 months. These features can be understood as the mutual 
interference between an assembly of discrete features in the scattered brightness distribution. 
By model-fitting to the observed arcs at one epoch we obtain a "snap-shot" estimate of the 
scattered brightness, which we show to be highly anisotropic (axial ratio > 10 : 1), to be centered 
significantly off axis and to have a small number of discrete maxima, which are coarser the speckle 
expected from a Kolmogorov spectrum of interstellar plasma density. The results suggest the 
effects of highly localized discrete scattering regions which subtend 0.1-1 mas, but can scatter (or 
refract) the radiation by angles that are five or more times larger. 

Subject headings: ISM: general — scattering — plasmas — pulsars: individual (B1133+16) 



1. Introduction 

The scattering of the radio waves as they prop- 
agate through the ISM has long been a clue to 
small scale inhomogeneities in the ionized inter- 
stellar plasma. The scattering is seen in the dy- 
namic spectrum of pulsar signals as random mod- 
ulations in frequency and time. Islands of inten- 
sity in the dynamic spectrum, often referred to 
as "scintles" , are sometimes modulated by a criss- 
cross pattern. By Fourier analyzing the dynamic 
spectrum, Stinebring et al. (2001) discovered the 
remarkable phenomenon of parabolic arcs that un- 
derlies the criss-cross substructure. In this pa- 
per, we analyze parabolic arcs from archival ob- 
servations of PSR B1133+16 originally reported 
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by Gupta, Rickett, & Lyne (1994). 

The two-dimensional Fourier spectrum of the 
primary dynamic spectrum is often referred to as 
the secondary spectrum. However, as described 
by Cordes et al. (2006), it can also be regarded 
as the "differential delay-Doppler" distribution, 
whose coordinates are time delay (/„) and Doppler 
frequency (or fringe rate ft), being the variables 
conjugate to radio frequency v and time t. The 
theory of the arcs is described in that paper, which 
henceforth we refer to as CRSC, and is also dis- 
cussed by Walker at al. (2004). 

The basic arc phenomenon is remarkably sim- 
ple, even though it lay undiscovered for over 30 
years. When scattered waves arrive at the ob- 
server from two angles they interfere and cause 
fringes in intensity. The observer is moving rela- 
tive to this interference pattern and sees a sinu- 
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soidal variation in both time and frequency, which 
thus appears as a pair of points in the secondary 
spectrum. The coordinates of these points are the 
differences in group delay and Doppler frequency 
between the two scattered waves. The Doppler 
shift depends on the relative velocity between the 
source and observer along the (scattered) lines of 
sight. Hence the difference in Doppler shift of 
waves scattered along two different paths through 
the ISM depends only on the angles of scatter- 
ing and the transverse velocity of the observer and 
pulsar relative to the ISM. It is the beating of these 
slightly different Doppler shifts that makes the 
fringe vary in time and so their difference equals 
the fringe rate ft- There is a quadratic relation 
of delay to fringe rate because, while the delay 
depends on the square of the angle of scattering, 
the Doppler shift depends linearly on the angles. 
When one of the waves is unscattered the result is 
a simple parabola f v oc ff. 

The observational properties of pulsar scintilla- 
tion arcs have been described by Hill et al. (2003), 
who observed arcs over a broad range of wave- 
lengths and found the curvature to be remark- 
ably constant and to follow the expected scaling 
as wavelength squared. Arc observations from five 
pulsars are shown by CRSC, who summarize the 
various arc phenomenon. The results show many 
deviations from the simple arc f v oc / t 2 . A partic- 
uarly striking form is the occasional appearance 
of "reverse arclets" which are displaced from the 
origin and whose curvature is equal in magnitude, 
but reversed in sign, to that of the main (forward) 
arc. Hill ct al. (2005) reported on a set of remark- 
able stable reverse arclets from PSR B0834+06. 

We present results for the evolution of the arcs 
seen in PSR B1133+16 at 408 MHz in archival 
observations spanning 16 months. On two of 
the days observed with short time resolution the 
dclay-Doppler spectrum showed several discrete 
arcs with reversed curvature, similar to those 
shown by Hill et al. (2005). We concentrate on 
fitting a model to one of those days. From an as- 
sumed brightness distribution we calculate the ex- 
pected delay-Doppler spectrum using the screen 
theory of CRSC. Then we iteratively adjust the 
parameters of the model to fit to the observed 
spectrum. Although the fitting is not unique we 
are able to explain the reversed arcs and obtain 
an estimate of the scattered brightness distribu- 



tion. This model distribution is then compared 
with that expected from a Kolmogorov scattering 
medium. We finish with a brief discussion of what 
structures in the interstellar medium might be re- 
sponsible. 

The technique we use, which is model-fitting 
to match the observed secondary spectrum, dif- 
fers from the inversion described by Walker and 
Stinebring (2005), who perform the match to the 
primary dynamic spectrum. We discuss the differ- 
ences between the two methods in £13.11 

2. Observations and Data Analysis 

The observations were made with the 76m tele- 
scope at Jodrell Bank between September 1984 
and January 1986, as described by Gupta, Rickett, 
& Lyne (1994). An auto-correlation spectrometer 
was centered at 408 MHz with 256 channels of 20 
kHz covering 5 MHz. The spectrometer was gated 
synchronously with the pulsar period and spectra 
were averaged in windows on and off the pulse. 
The averaging time was normally 30 seconds and 
on two occasions it was shortened to 10 seconds. 
The duration (T b s ) of most observations was from 
one to four hours. The on and off pulse spec- 
tra were corrected for the one-bit auto-correlation 
normalized and subtracted to create a primary dy- 
namic spectrum of pulse intensity S(v,t). 

The left hand panels of Figure Q] shows S(i>, t) 
for 8 and 9 September 1984 (which we call days 
and 1 of the observing sequence); they exhibit 
crisscross structure with a predominance of scin- 
tles that drift toward higher frequency with time. 
With the 10 second sampling one can see fine scale 
fringes drifting in the same direction. 

The right panels are the corresponding sec- 
ondary spectrum - S2(fv,ft) or delay-Doppler 
spectrum, which is the squared magnitude of the 
Fourier Transform of the primary spectrum: 

S(y,t) — > St(/„,A) (1) 
S 2 (fuJt) = \S\U,f t )\ 2 . (2) 

Here S^(f v ,ft) is the two dimensional discrete 
Fourier transform of the sampled primary spec- 
trum S(v,t). With time resolution of 10 seconds 
and frequency resolution of 20kHz, the Nyquist 
frequencies are /t,Nyq = 3 cycles min -1 (cpm) and 
/y,Nyq = 25.6 cycles/MHz (or 25.6 //sec). The as- 
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Frequency v (MHz) Delay f (microsec) 



Fig. 1— Top Left: Dynamic Spectrum of PSR B1133+16 on 8th September 1984 (day 0) with sampling 
time of 10 seconds. Top Right: Secondary Spectrum of the same data. The plot is limited to ±2 cpm in 
ft, since no signal is visible in the range 2-3 cpm. The grey scale is logarithmic from the mean noise level 
to black at a level 30dB higher, as shown by the tablet to the right. The plot shows a set of reversed arcs 
grouped together forming a forward parabolic arc predominantly visible for negative ft- Bottom Left and 
Right: Same displays for 9th September 1984. Note that the different time scale from that above; both show 
narrow fringes that drift in the same general sense. Both secondary spectra are similar, but the details of 
the reverse arcs have changed over one day. 



sociated spectral resolutions are 5ft = 1/T ^ S and 
&fv = 0.2^isec. 

The S2 plots show structure which is strongest 
near the origin and extends out to large delays 
(f u ) in the lower part of the panel (ie for negative 
ft). We call this the forward arc, since it is curved 
toward the /„ axis. It consists of discrete curves, 
which we call reverse arcs, since they are curved in 
the opposite sense. The apexes of the reverse arcs 
lie approximately along the parabola f v oc ff . On 
close inspection one can see several nearly parallel 
forward arcs. The top and bottom panels of Figure 
[1] are separated by 27 hrs and show very similar 



behavior in general, but the individual reverse arcs 
have clearly changed substantially. The main goal 
of our paper is to investigate these reverse arcs 
and to find a model for the interstellar angular 
scattering function that could explain them. 

2.1. Time Evolution of the Arcs 

After rapidly sampled observations were made 
on September 8th and 9th 1984, a slower stan- 
dard sampling interval of 30 seconds was adopted 
for monitoring the ISS of PSR B1133+16. Figure 
[5] shows S*2 plots from twelve observations span- 
ning 481 days. Days 0-31 show the same general 
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Fig. 2. — Selected Secondary Spectra of PSR B1133+16 from September 1984 to January 1986 in a loga- 
rithmic grey scale over a 30 dB dynamic range, as shown by the tablet in Figure [1] 



character of a V-shaped "valley" with little power 
along the delay axis, flanked by diffuse power often 
appearing as a set of discrete reverse arcs, which 
can be asymmetric in f t . Then day 45 shows a 
simpler forward arc which fades remarkably over 
days 96 and 142. A strong forward arc appears on 
day 194 and by day 288 multiple reverse arcs sim- 
ilar to days 0-31 have returned. Day 415 is similar 
to day 194, but day 481 has weak and indistinct 
arc structure. These remarkable changes in the 
appearance of the secondary spectra are evidence 
that the line of sight crossed through particular 
localized structures, as we discuss in § 4. 

With 30 second sampling /t,Nyq — 1 cpm, and 
so features beyond ±1 cpm can be aliassed into 
the figure. This is evident by comparing the plots 
for days and 1 in figures ((2|) and {]]). For exam- 
ple on day 1 figure ([2|) shows a linear feature at 
ft ~ 0.75 ± 0.15 cpm and f v ~ 18^sec. This is an 
alias from f t ~ 0.75 — 2.0 = —1.25 cpm, visible in 
the lower right panel of Figure [1] at /„ ~ 18/xsec 
with its apex near f t = —1.4 cpm. The figure ([1]) 
data were a 23 minute observation immediately 
following the 30 min observation figure ([2|). This 
discrete reverse arc changed very little in 40 min- 
utes. Though it clearly changed between day 
and 1. Days 9 and 16 show fuzzy reverse arcs for 
both positive and negative ft, which are probably 
confused by aliasing, but give evidence for an evo- 
lution in the asymmetry as well as in the reverse 
arc features. 

In one of the observations of PSR B1133+16 
shown by CRSC the secondary spectrum has two 
very well defined discrete forward arcs. Putney 
and Stinebring (2006) have found multiple arcs 
in six pulsars. For B1133+16 they found the 
arc curvatures to fall near four discrete values, 
but that the arc strengths were quite variable. 
They suggest that these are due to variable am- 
plitude scattering at four discrete distances from 
the Earth. Consequently, we have estimated the 
curvature a (defined by f v = a/ t 2 ) from each of 
our observations scaled to 1 GHz. The mean of 
our values over 1984 September to 1986 January 
was a = 6.7 x 10 -3 sec 3 with an rms scatter of 
0.6 x 10 -3 sec 3 , which is consistent with one of 
their values at (6.3 ± 0.1) x 10 -3 sec 3 . We do not 
detect multiple forward arcs, but note the signal- 
to-noise ratio in our data is lower than in most of 
Stinebring's data. We also note that the curvature 



estimation becomes more difficult in the presence 
of reverse arclets, and in such cases we trace the 
arc by the location of their apexes. 

3. Modeling of arcs from PSR B1133+16 
3.1. Theory 

We use the theory of arcs in the secondary spec- 
trum as given by CRSC. Section 3.3 of that pa- 
per analyzes S2 in terms of a "brightness func- 
tion" B(9 x ,9 y ) that represents the instantaneous 
strength of contributions to the intensity at a given 
observing point as a function of the angle of arrival 
9 = (6 X) 6y). The theory provides insight into the 
basic interference effects that give rise to the arc 
phenomenon, but it does not include the disper- 
sive delays introduced by propagation through the 
irregular plasma density in the ISM. In fl3.5l we use 
a full electromagnetic simulation to include such 
effects. 

Consider a pulsar at distance D whose radia- 
tion is scattered by a screen at a distance sD from 
the pulsar. The observer sees a superposition of 
interference fringes from each pair of apparent an- 
gles of arrival Q\ and 82- The fringes contribute 
to a single pixel in S^C/v, ft) if the angles obey: 




where c is the speed of light and A is the central 
wavelength observed. Vj SS is the effective trans- 
verse velocity relative to the scattering region of 
the line of sight from pulsar to observer, which 
both may be moving (Cordes and Rickett, 1998); 
without loss of generality we have chosen the ori- 
entation of Vi ss to define the direction of X ■ 

Equation (8) of CRSC expresses S<z{f u , f t ) as a 
four-fold integral over all pairs of angles, subject 
to the constraints of equations [3] and |U The con- 
straints allow the result to be reduced to a two- 
fold integral as given by CRSC equations (9) or 
(B4 & B5). We used the latter, since it is more 
convenient numerically, and we rewrite it here in 
a slightly different form as s 2 {q 7 p) which is the 
delay-Doppler spectrum in terms of normalized de- 
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lay and fringe rate: 
2cs 



P = fv 



D(l- S )6 2 ' 



q = ft JTT ( 5 ) 



9q is a normalizing angle, typically chosen to be 
the angular radius of the diffractive brightness dis- 
tribution. 

The result is 



s 2 (p,q) 



OO pQC 
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q j- 

p , q i 2 



(0) 



t± = —±--—{V2-Vi) 
^ 2q 2 2q Ky y ' 



where b(x 1 y) is the brightness function in terms 
of normalized angular coordinates x = |£ , y = -j- , 
and the constraints in equations ([3] & 2| become: 



x 2 



P 
x 1 . 



(x 2 



2 x 1 2 ) + (y 2 2 -y 1 2 ) 



(7) 



There is also a second mirrored relation obtained 
by mapping p => — p, q =>• — g. 

Hence we can describe an arc as the locus of 
points in p, q due to interference of a narrow dis- 
crete component in b(x,y) at x p ,y p with the re- 
maining extended continuum in b, as is done in sec- 
tion 3.4.3 of CRSC. A primary forward arc p > q 2 
comes from an undeviated discrete component in 
b interfering with the extended continuum in b. 
This is the condition under weak scintillation de- 
scribed by CRSC. The boundary of this primary 
arc is f v = af 2 (or p = q 2 ). In turn a is related to 
the physical parameters via equations H] and [3J 

In general we substitute x p ,y p for 22,2/2 in 
equation ((7|) to obtain a relation between p and 
q depending on X\ and yi, which vary over the 
continuum in b. The special case where b is nar- 
row in one dimension and extended along a line 
yi = 2itan0, at an angle <j) to the i-axis, al- 
lows us to eliminate both x\ and y\ to obtain the 
quadratic relation between p and q: 



P 



V 2 P -sec 2 (</>) {q-x p f 



(8) 



This reverse arc is a parabola with an apex at 
Pa = x 2 + y 2 ,q a — x p and negative curvature 
— sec 2 ((f)), which has unit magnitude if <fr = (i.e. 
scatter broadening parallel to the velocity). The 
corresponding forward arc with its apex at nega- 
tive p is obtained by reversing the signs of p and 



q in equation ([8]). Both sets of arcs are shown in 
Figure 4 of CRSC. 

Of course, in general the scatter broaden- 
ing function is not one-dimensional, as assumed 
above. Further, we have not discussed the power 
density of s 2 along an arc. This depends on the 
relative strength of the discrete component sup- 
posed to exist in b and on the shape and intensity 
of its continuum, as described by the integral ^ . 

3.2. Model for September 8th 1984 (day 
0) 

The parallax and proper motion of PSR Bl 133+16 
were measured by Brisken et al. (2002); they 
found its distance D = 0.35 ± 0.02 kpc and its 
transverse velocity V pm = 631 ± 36 km sec -1 . 
In general the scintillation velocity depends on 
the transverse velocities of the pulsar, the Earth 
and the scattering medium, but since our pulsar 
is moving much faster than the Earth and the 
medium, we have V; ss = V r pm (l — s). Hence the arc 
curvature, as given by equation (8) of Stinebring 
ct al. (2001), becomes: 



DX 2 



2cV 2 m 1 - s 



0.0244- 



1 - s 



(9) 



at our observing wavelength. The arc curvature 
a = 6.7 x 10 -3 sec 3 , given in section 2.1, is scaled 
to 1 GHz. When scaled back to our observed fre- 
quency of 408 MHz, we have a = 0.040 ± 0.003 
sec 3 . This results in s = 0.62, giving a screen at 
133 ± 7 pc from the Earth. Putting the mean 
value into equations [3] and [4] gives the mappings: 



t.cpm 



= 0.151 



s) 



0.26 



□2 
2mas 



9 2 \ 
lmas / 



(10) 



The secondary spectrum from September 8th 
1984 shown in the upper right panel of Figure Q] 
exhibits a series of reverse arcs extending to de- 
lays of 25 /^sec, which are visible predominantly 
at negative / t . We choose these data to model 
quantitatively, because of these unusual asymmet- 
ric discrete reverse arcs. Consider, for example, 
the faint reverse arc with its apex at f v = 14 [isec 
and ft = —1.1 cpm, which is visible over the ap- 
proximate range —1.6 < ft < —0.8 cpm. 

A discrete reverse arc is formed by interference 
of a point component in the scattered brightness 
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function with an anisotropic extended component. 
To simplify the analysis consider the normalized 
brightness b(x, y) to be narrow in y and extended 
in x, i.e. parallel to the ISS velocity (i.e. = 
in equation©). Thus we write b(x,y) = [a p 8(x — 
x p ) + b e (x)]S(y), where b e represents the extended 
component and we use delta functions to model 
the narrow structures. The interference between 
waves from x p and x gives q — x p — x. Thus, if 
b e {x) extends from x_ to x + , the reverse arc will 
have its apex at p a = x 2 , q a = x p and be visible 
over X- — x p < q < x + — x p . 

In this model y — and the apex lies on the ba- 
sic arc. Hence substituting f v and ft from above 
into equation © gives s = 0.63, which is consis- 
tent with the value 0.63 obtained above. Using 
the mapping in equation 1 101 the observed range 
~ 1-6 < ft _ —0.8 cpm for the reverse arc, corre- 
sponds to interference of a discrete component at 
-7.5 mas interfering with an extended component 
from -2.7 mas to 2.7 mas. In a similar way for 
each discrete reverse arc on day we can find a 
corresponding pair of components in the scattered 
brightness (a point and extended component along 

e x ). 

Following this idea we model the brightness dis- 
tribution (in normalized coordinates) as a sum of 
Gaussian functions: 

b(x, y) = ^T Ait'^ 1 ^ (11) 

i=l 

where Ai is the intensity, a X i describes the width 
of x, Uyi describes the width of y, x i is the offset 
in x, and y Q i is the offset in y. We made an ini- 
tial model in the manner described above for the 
most obvious reverse arcs in Figure [TJ As more 
components are added the extended component 
eventually becomes the sum of the discrete com- 
ponents. Then we tried to optimize the model by 
minimizing the sum of the squared errors between 
the observed S2 and the model obtained by a nu- 
merical integral of equation 

3.3. Fitting Process 

We chose to work in normalized delay p, fringe 
rate q and angle x, y and used a normalizing an- 
gle 8 = 3.1 mas, with the constants evaluated as 
discussed above. The fitting was done over two 
regions where the arcs are visible in the top right 



panel of Figure [T] a lower rectangle -2-0 cpm by 
0-25 fj,s and a smaller upper rectangle 0-1 cpm by 
0-5 fjs. 

We estimated the noise background in from 
the average noise level in £2 over the regions not 
included in the fit, and this constant was added to 
the model. We defined x 2 in a standard fashion: 

^,2 _ [S^ifv^ii ft,j) ~ ^2,mod,i,j] 2 

a i,j 

where a it j = S2,mod,i,j- At eacn P oint £2 (/*,», f t ,j) 
is an estimate of the secondary spectrum, which 
has a standard deviation cr, j equal to its mean 
value (following exponential statistics). Once the 
model is close to an optimum value it provides an 
estimate of <7jj. Thus we used the model as an 
estimate for aij rather than the observed value. 

The fitting was done in a semi-automated fash- 
ion. Starting with separate optimization of the 
parameters describing each narrow component re- 
sponsible for the more obvious discrete arcs, we 
minimized \ 2 over a rectangle around each apex. 
Table Q] lists the coordinates of the apexes of six of 
the most obvious reverse arcs. With five parame- 
ters for each Gaussian the number of parameters 
to fit grows quickly as components are added. To 
simplify we started with all components centered 
on the x-axis (ie y Q i — 0). With the position x oi 
of a component fixed we explored the influence of 
<j X i and a y i . For the arcs at large delays we found 
that a x i had to be quite narrow to agree with the 
observed arc thickness in delay, but that <r y i could 
be substantially wider. As a result we did not op- 
timize for the y-position or width, but set them all 
on the x axis with a y i — 0.3 mas. This left three 
parameters to be optimized for each component 
(amplitude, s-position and width). However, as 
discussed above the extent in q is governed by the 
x- position and width of the associated extended 
component, thus they too had to be adjusted. 

Close inspection of the observed S2 shows nar- 
row parallel forward arcs crossing the reverse arcs. 
As can be seen from Figure 4 of CRSC, these are 
extensions of the reverse arc to positive f t values 
which are reflected about the origin. Thus they 
show interference of a discrete component (x p < 0) 
with components of the scattered brightness ex- 
tended over a wide range, extending to x values 
more negative than x p . This emphasizes that in 
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fact the entire set of arcs is caused by the mutual 
interference of all possible pairs of points in b(x, y). 

As noted above the observations show parallel 
forward arcs (for negative ft), but they also in- 
clude some that appear as fainter (whiter) stripes. 
Such fainter stripes come from a discrete min- 
imum in b(x,y), whereas the brighter arcs are 
from discrete maxima. We modelled this effect by 
putting dips into the brightness function, which 
we achieved by having two offset broad extended 
components with a minimum in b between them. 
These minima improved the fit slightly over a 
model without the minima. Thus in choosing how 
many and where to place new components, we 
were guided by the visual appearance of the ob- 
served S2 in comparison with the model. 

Table [2] gives a list of the parameters for the 
9-component model for b(x,y) (scaled into milli- 
arcseconds). Components 1-6 are relatively nar- 
row [a x < 0.3 mas) and components 7-9 are rel- 
atively wide. It is the interference of the narrow 
with the wide components that is responsible for 
the 6 most obvious discrete reverse arcs. Their 
position and width in Doppler frequency governs 
the range in x of the brightness function. 

We optimized the model to minimize \ 2 with 
regard to the three variable parameters of the 9 
components. The model can be seen in the upper 
panel of Figure El to be compared with the upper 
right panel of Figure [TJ However, the value of \ 2 
was about five times the number of points in the 
fit, indicating that the model is far from optimum 
even though it was a local minimum. The largest 
discrepancy was traced to the region near the ft 
axis (from -1 to -2 cpm), in which the observed 
power level is elevated but the model is not. The 
higher level in the observations appears to be due 
to effects from broadband variation of the pulse 
strength (intrinsic pulse modulations or impulsive 
interference). Consequently, we estimated this ex- 
tra power by averaging S2 in (from -2 to -3 cpm) 
as a function of delay, and added it to the model. 

This refinement caused little change in the best 
parameters but reduced the minimum x 2 to about 
twice the number of data points. Visually, the 
model resembles the observations in the asymme- 
try versus ft, the presence of 6 reverse arcs, the 
location of their apexes, and their width in de- 
lay (/„). Given the success of the fit from a visual 
perspective, we decided not to attempt a better fit 



by introducing extra components, since we would 
just be refining the model for a particular snap- 
shot of the scattering function b. Further, since 
our method of fitting relies on the visual impres- 
sion in choosing how many components and their 
initial parameters, the result is not unique. 

The upper panel of Figure [3] shows the model 
secondary spectrum. It corresponds to the bright- 
ness model in Table [5] with one modification. In 
order to represent the "speckle" effect in a snap- 
shot of the scattered brightness we multiplied the 
sum of Gaussians by a random function of x (uni- 
form distribution), before computing the integral 
in equation (J6|). The effect is to cause a set of 
fine nearly parallel forward arcs, somewhat remi- 
niscent of the observed fine structure in 5*2. The 
model brightness is plotted logarithmically in the 
upper panel of Figure 01 Note the asymmetry and 
anisotropy in its distribution of discrete compo- 
nents along the axis parallel to Vi SS . 

The effect of rotating the brightness by an an- 
gle to the cc-axis increases the magnitude of the 
curvature in p, q plane, as can be seen from equa- 
tion JHJ. However, this could be compensated by 
a change in the screen distance parameter s. Thus 
we cannot determine both s and <f> and leave our 
results displayed with = 0. 

We now consider what constraint can be placed 
on the y-extent of the brightness function (b(x, y)). 
When y p and y\ are non-zero, we obtain the a 
version of equation © in which y^ is replaced by 
yp — y\. Let the extended and point components 
be centered on the x-axis with widths a y and a yv , 
respectively. Then these ranges broaden the thick- 
ness of the reverse arc in delay (ie p is broadened 
by the larger of and cr^ p . The widths of the re- 
verse arcs in the upper right panel of Figure [1] are 
about ±0.05^sec. With the scalings from equa- 
tion (fT0| . this constrains the rms y- widths to be 
less than about 0.44 mas. We set them (somewhat 
arbitrarily) to be 0.3 mas. 

3.4. Discussion of the Model 

Although the parameters we used gave us arc 
and subarc structures that resemble the particular 
PSR Bl 133+16 data, we do not imply that it is the 
only set nor is it the best set. We have produced 
one model for the brightness distribution that re- 
produces the visual form of the observations and 
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is optimized within a limited range of the param- 
eters. 

From this we conclude that the asymmetri- 
cal form of the observed secondary spectrum and 
the presence of narrow discrete reverse arcs can 
be successfully explained by an asymmetrical and 
anisotropic brightness function with a number of 
narrow discrete features. The constraint on the 
widths in 9 y is important since, with an x -width 
of about ±5 mas, it constrains the effective ax- 
ial ratio of the scattered brightness function to be 
greater than about 10. However, an interesting 
feature of the model is that the individual features 
can be wider in 9 y than in 9 X ; in other words they 
may be anisotropic in the orthogonal direction. 

As shown by CRSC, an absence of power along 
the delay axis in the secondary spectrum is caused 
by anisotropic broadening, enhanced along the di- 
rection of the scintillation velocity. All of 12 of 
the secondary spectra plotted in Figure [5] show 
no power along the delay axis. Hence the same 
basic anisotropy in the scattering persisted over 
about 500 days, but the detailed form of the scat- 
tered brightness changed substantially over days 
and weeks. 

3.5. Screen Simulation 

CRSC used the simulation code of Coles et 
al. (1995) to simulate secondary spectra of waves 
scattered by a dispersive phase changing screen. 
In normal usage the method synthesizes one real- 
ization of a stationary random phase screen, which 
modulates a plane wave. The emerging complex 
field is then propagated a distance z to the ob- 
server, by transforming into wavenumber and mul- 
tiplying by the appropriate (Fresnel) filter and 
back transforming to the spatial domain. The in- 
tensity is saved along a straight line track through 
the two-dimensional complex wave field. The re- 
sult simulates a time series as an observer moves 
relative to a scattered wave field. A dynamic spec- 
trum is a set of such time series, created as the ra- 
dio frequency of the incident wave is incremented. 
The secondary spectrum is computed in the same 
way as for the data. 

CRSC simulated the screen as a single two- 
dimensional realization of a stationary random 
process following a given wavenumber spectrum. 
They showed results for a Kolmogorov screen. 



When anisotropy was included, some results had 
multiple discrete forward and reverse arcs. How- 
ever, the results differed from our observations of 
PSR B1133+16 in that the discrete arcs were con- 
siderably finer and distributed symmetrically in 
fringe rate (ft)- The addition of a linear gradient 
in the phase across the screen (i.e. a simple refrac- 
tion) caused the secondary spectrum to be asym- 
metric in fringe rate. They showed that the phase 
gradient needed to be large enough to shift the 
scattered image by at least its own radius. This is 
a much larger gradient than is expected randomly 
for a medium with a Kolmogorov spectrum. 




1 5 10 , 15 20 25 

t v (microsecond) 



Fig. 3. — Top: Secondary spectrum computed 
from the brightness distribution of the upper panel 
of Figure 31 by a numerical integration of equation 
(J5]) ■ Grey scale is logarithmic with dynamic range 
of 36 dB. The fit was done over the two rectangu- 
lar regions shown; outside those regions the noise 
from the data is plotted. Bottom: Secondary spec- 
trum from the screen simulation for the same day 
(see text) with dynamic range of 30 dB. 

We modified the screen generation code to make 
the scattered brightness like that in our model 
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(i.e. the upper panel of Figure [4]). This modifi- 
cation required an iteratative process, as follows. 
The field emerging from the phase screen at posi- 
tion s has unit amplitude and is phase modulated 
F(s) = exp[j0(s)]. Its two-dimensionsal Fourier 
transform F^(k) is the spectrum of the field ver- 
sus wavenumber n, which maps simply to the scat- 
tered angle as k — k6 with k = 2tt/X. Our goal is 
to find a screen phase <fi(s), for which the squared 
magnitude |Ft| 2 follows our model for the scat- 
tered brightness B(0). We started by specifying 
the mag nitude = ^B{6/k) and assigned it 
random phase. Then back transformed to give 
F(s), which had both amplitude and phase mod- 
ulation. A revised F(s) was obtained by normal- 
izing each point to unit magnitude. That revised 
F was re-transformed and the resulting F^ was 
normalized to make its magnitude y/B(0 /k) and 
keep its phase. When this cycle was repeated the 
amplitude modulation in F(s) is decreased. We 
iterated until the residual amplitude modulation 
no longer decreased (~ 100 cycles). The phase of 
this final F(s) was then extracted as the screen 
phase for the normal simulation. This method 
is based on the concepts developed by Gerchberg 
and Saxton (1972). As in the normal simulation 
this screen phase was scaled proportional to radio 
wavelength, to represent the plasma dispersion. 

The lower panel of Figure [3] gives the resulting 
secondary spectrum. It has the same basic fea- 
tures as the simple model in the upper panel, but 
shows more pronounced substructure in the form 
of fine arc modulations. The importance of this 
result is that the basic arc forms are unchanged 
by the inclusion of dispersion in the screen phase, 
whereas dispersion is not included in the physics of 
equation ([6]). A further result is that we have ob- 
tained a realization of the phase screen itself. This 
is then an estimate of the phase imposed by the 
ISM during the observations on our day 0, and we 
can ask what it tells about the density structure 
in the ISM. 

3.6. Comparison with Reverse Arcs from 
PSR B0834+08 

Hill, et al. (2005) reported how four discrete 
reverse arcs observed from PSR B0834+16 moved 
along the main forward parabola over 25 days. 
The Doppler frequency of the apex of each arc 
followed a linear track against time; furthermore 



the slopes of all four tracks were equal to the slope 
predicted by assuming that each reverse arc was 
scattered from a fixed structure in space scanned 
by the known pulsar proper motion. Their results 
suggest that the reverse arcs appearing on succes- 
sive days are not unrelated, but are in fact due to 
the same physical structures. 

We have attempted to test this idea with the 
reverse arcs observed on day and 1 of our data 
(PSR B1133+16), which are 27 hours apart. For 
the 6 most obvious reverse arcs on day we calcu- 
lated the shift expected due to the pulsar proper 
motion in 27 hours and compared the new loca- 
tions to those observed on day 1. The expected 
angular shift is S6 X ~ 5tV pm /D, which we sub- 
stitute into equation [T0l and tabulate the ft and 
f v positions before and after this proper motion 
shift (see Table |3]) . Comparing the shifted coor- 
dinates with those observed on day 1, we found 
some agreement. For example, the easily visible 
reverse arc no. 4 from day should should move 
to ft — -1.25 cpm and /„ = 18.75 fj,s which is close 
to a prominent arc on day 1 (lower right panel 
of Figure Q}. The comparison is complicated by 
the fact that day 1 shows a higher density of re- 
verse arcs making it hard to recognize particular 
structures. We conclude that these data are not 
of sufficiently high signal to noise ratio to confirm 
or reject the hypothesis that the arc positions are 
simply shifted by the pulsar proper motion. 

3.7. Relationship to W-S Inversion Algo- 
rithm 

We now compare our model fitting with the 
method of Walker and Stinebring (2005) (WS) 
who estimate the scattered brightness function us- 
ing an fitting algorithm to the dynamic spectrum. 
This has great appeal since it avoids the messy 
task of deciding what components of the secondary 
spectrum to include in a model fit. Their method 
represents the full electric field as a sum of "scat- 
tered waves" , which are successively identified and 
subtracted in order to match the observed primary 
dynamic spectrum to the model. It is somewhat 
analogous to the CLEAN algorithm in image re- 
construction. 

WS consider the scattered electric field U{v,t) 
at frequency v and time t as a sum of basis func- 
tions exp[i2ir(vft + f„t)], which they call scattered 
waves, though the formulation is not derived from 
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scattering theory. We note that v and t are relative 
to the center of the observing band and integra- 
tion time and that WS use the symbols r for /„ 
and to for f t . 

The dynamic spectrum S(v, t) is then modelled 
by \U(v, t)\ 2 and WS designed a procedure to es- 
timate the phasor (U(f v ,ft)) for each scattered 
wave to best match the model with the observed 
dynamic spectrum. The results are given as both 
amplitude and phase for each differential delay f v 
and Doppler f t . They applied the method and 
displayed the amplitudes for a particular 327 MHz 
observation in the sequence from PSR B0834+06 
(Hill et al., 2005), in which there were several well- 
defined reverse arclets. They discuss possible ap- 
plications of the method but do not relate the re- 
sults to the scatttered brightness B(9 x ,0 y ). 

However, Asplund et al. (2004) used these same 
results and mapped them to the scatttered bright- 
ness by assuming that each component U(f v ,ft) 
was the result of interference between an undevi- 
ated (reference) wave and a wave scattered at an- 
gle (9 X , 6 y ). This assumption is valid under condi- 
tions of weak scattering, as shown by CRSC, but 
it does not seem justified for the 327 MHz ob- 
servations of PSR B0834+06. As WS note, the 
effective reference wave in their algorithm starts 
out as undeviated, but is modified as each scat- 
tered wave component is identified. Nevertheless 
the summing of the "scattered waves" does include 
the mutual interference between the components 
of the basis functions as well as their interference 
with the reference component. 

One can show that the secondary spectrum is 
the square of the auto-correlation of U{f v ,ft), a s 
follows: 

S a (U, ft) = I EE U[F v ,F t ) U*{F V + U, F t + f t )\ 2 (13) 

Nevertheless, it seems that more work is needed on 
the theoretical basis of the algorithm to consider 
the circumstances under which the weak scintil- 
lation mapping into the angle domain is properly 
justified. 

While the method we have used is cumbersome 
to implement, it is more straightforward in that 
our computed model directly sums the waves in 
the domain of (8 x ,9 y ) including their mutual in- 
terference. 



4. The Density structure in the Interstel- 
lar Plasma 

Most interpretations of radio scintillation have 
been formulated as constraints on the wavenum- 
ber spectrum of the plasma density. The Kol- 
mogorov spectrum has become a default model, 
even though several observations point to more 
fluctuation power on wavelengths longer than the 
size of the typical radio scattering disc (a few AU) . 

The general concept proposed to explain the 
discrete reversed arcs is that substructure in the 
scattered brightness b, analogous to speckle in a 
scattered image, causes substructure in the sec- 
ondary spectrum, that takes on a form of points 
lying on intersecting reverse and forward parabolic 
arcs, traced by equation (JSJ) and its reflection in 
the origin, with intensity at each point governed 
by the product of the brightness of the two in- 
terfering components. The plots in Figure 3] are 
estimates of a single realization of the scattered 
brightness with a particular sub-structure charac- 
teristic of the observing conditions. We now con- 
sider what these images tell us about the density 
structure in the interstellar plasma. 



-10 -5 5 

K (mas) 




« (mas) 



Fig. 4. — Top: Logarithmic plot of the model 
brightness distribution obtained by fitting to ob- 
servation on day in Figure [T] (dynamic range 
50 dB). Bottom: brightness distribution obtained 
from the full screen simulation (see text). 

The first point is that the scattered image is 
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highly elongated with an axial ratio of 10:1 or 
more. Further, the axis of elongation is roughly 
parallel to the proper motion direction of the pul- 
sar (within, say, ±20°). This means that the spa- 
tial density structures are elongated in a direction 
roughly perpendicular to this direction. 

The second point is that the images are not 
centered on the origin. We see two possible expla- 
nations: cither refraction by a large scale trans- 
verse gradient in the column density of electrons; 
or highly localized scattering caused by a region 
that is centered to one side of the line towards the 
pulsar. 

Refraction seems likely since similar refraction 
has been invoked to explain the persistent drift- 
ing in frequency of scintillation features in the dy- 
namic spectra of some pulsars (e.g. Shishov, 1974; 
Hewish, Wolszczan and Graham, 1985; Gupta, 
Rickctt and Lyne, 1994). The basic idea is that 
a refracting region deflects the scattered (and un- 
scattered) waves which causes a lateral shift in the 
scintillation pattern. Since in a plasma the angle 
of refraction is strongly dependent on the observ- 
ing frequency, the frequency for a given maximum 
of intensity will drift in time. Such a drift will be 
seen as an asymmetry versus Dopplcr frequency in 
the power distribution of the secondary spectrum, 
this effect was confirmed by simulations reported 
by CRSC. 

The alternative explanation is that there are 
highly localized regions of very irregular (turbu- 
lent?) plasma, that can deflect (ie scatter or re- 
fract) the radio waves into our line of sight. This 
implies that thay must be small enough to sub- 
tend an angle at the observer smaller than their 
maximum angle of deflection. The angles of ar- 
rival of these waves will be determined by how far 
the enhanced turbulence is from the direct line of 
sight and they will not depend systematically on 
frequency. 

A third point about the scattered images is 
that they are "blotchy", but the blotches have a 
lower filling factor than would be expected for the 
usual speckle in an angular spectrum, (see the 
simulations shown by CRSC, in which there is a 
multitude of reverse arcs due to speckle with a 
high filling factor from a Kolmogorov medium). 
We conclude that the blotchy nature of our scat- 
tered image is not compatible with a uniform Kol- 
mogorov random medium. Our results do not, on 



their own, allow us to determine the nature of the 
medium. However, the fascinating measurements 
by Hill et al. (2005) for pulsar PSR B0834+06 
provide an important constraint. They found that 
the frequency-scaling of the positions of the re- 
verse arcs was not compatible with frequency de- 
pendent refraction and implies the alternative of 
highly localized regions of enhanced scattering as 
the cause of each reverse arc. Under this scenario 
the regions have a linear size given by their ap- 
parent angular size times the scattering distance 
- D(l — s) ~ 140 pc. With angular extents of 
0.1-1 mas their maximum implied linear sizes are 
0.014-0.14 AU. 

Considering that the discrete reverse arcs were 
only prominent in some of the observations (ie 
days 0-31 and 288), the regions of localized en- 
hanced scattering are irregularly distributed in the 
ionized ISM. In 31 days the pulsar moves about 10 
AU transverse to the line of sight. Thus the persis- 
tence of the reverse arcs through 31 days suggests 
that the localized regions of enhanced scattering 
are clumped over about 10 AU. 

If localized scattering regions are indeed respon- 
sible for the discrete reverse arcs, the next ques- 
tion is what are these regions physically. However, 
this is beyond the reach of our present investiga- 
tion, and we refer interested readers to Rickett, 
Lyne and Gupta (1997) and the proceedings of the 
workshop on "Small Ionized and Neutral Structure 
in the ISM", held in Socorro May 2006, particu- 
larly the papers by Stinebring and Rickett. 

5. Conclusion 

The ISS of pulsar Bl 133+16 showed pro- 
nounced changes in the character of its secondary 
spectra in 12 days of observation spread over 15 
months. While these data which were centered at 
408 MHz were recorded over 20 years ago, they 
often show the parabolic arc phenomenon that 
was discovered quire recently by Stinebring et al. 
(2001). However, the basic parabolic arc varied 
substantially in its extent and shape. Discrete re- 
verse arcs were visible for about two months early 
in the sequence and then again after about nine 
months. There was a two month interval when 
the arc was barely detectable, corresponding to 
when the primary dynamic spectra had a very 
wide decorrclation bandwidth. 
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We concentrated our analysis on two days with 
clearly visible discrete reverse arcs, which were for- 
tuitously sampled rapidly enough to resolve them 
in Dopplcr frequency. Using the theory of CRSC, 
which ignores dispersion, we modelled these fea- 
tures and hence were able to estimate the scat- 
tered brightness distribution. We tested the ap- 
plicability of the theory, by simulating and elec- 
tromagnetic wave propagation code that includes 
dispersion, and found that the reverse arcs are not 
altered substantially by including dispersion. 

Our estimate of the scattered brightness func- 
tion is highly clcongatcd with an axial ratio of 
more than 10:1. It is asymmetrical about the ori- 
gin and has a small number of discrete maxima. 
These properties are not consistent with expecta- 
tions from a scattering region that behaves as a 
statistically homogeneous random medium with a 
Kolmogorov spectrum. When combined witth the 
conclusions from Hill et al. (2005), our results sug- 
gest that there are highly localized discrete scat- 
tering regions which subtend 0.1-1 mas, but can 
scatter (or refract) the radiation by angles that 
are five or more times larger. The asymmetry in 
the scattered image may well be caused by the 
chance configuration of these regions. 

We thank Andrew Lync for completing the se- 
ries of observations at Jodrell Bank in 1984-85 and 
we thank Yashwant Gupta for work on the origi- 
nal data reduction. We thank the National Science 
Foundation for support under grants AST 9988398 
and AST 0507713. 
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This 2-column preprint was prepared with the AAS IAT^X 
macros v5.2. 



Subarc Coordinates (Observation vs. Model) 



Observation Model 
Subarc # ft (cpm) f u ((is) ft (cpm) f v (/is) 



1 


-0.31 


1.2 


-0.32 


1.22 


2 


-0.58 


3.8 


-0.56 


3.82 


3 


-0.69 


5.4 


-0.67 


5.43 


4 


-1.09 


14.0 


-1.10 


14.22 


5 


-1.31 


20.0 


-1.30 


20.36 


6 


-1.44 


23.6 


-1.41 


23.86 



Table 1: The coordinates of the observed individ- 
ual subarcs (day 0, Figure 1: Top Right) and the 
modeled subarcs (Figure 3: Top). 
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Parameters for Model Brightness Distribution 



Component 



(i) 




[%oi ) 




(Voi) 


((Tyi) 


1 


0.418 


-2.18 


0.138 





0.32 


2 


0.273 


-3.85 


0.104 





0.32 


3 


0.127 


-4.59 


0.126 





0.32 


4 


0.327 


-7.42 


0.041 





0.32 


5 


0.2 


-8.88 


0.261 





0.32 


6 


0.236 


-9.62 


0.038 





0.32 


7 


1.0 


0.82 


1.04 





0.32 


8 


0.545 


-2.68 


1.418 





0.32 


9 


0.027 


-7.25 


2.111 





0.32 



Table 2: Components of the brightness distribu- 
tion defined by equation [Tl] where the angular 
units are mas. These values were used in Figured] 
to model the observations of day 



Subarc Coordinates (Day vs. Proper Motion) 



Before After 
Subarc # f t (cpm) f v (ps) f t (cpm) f v (/xs) 



1 


-0.32 


1.22 


-0.48 


2.78 


2 


-0.56 


3.82 


-0.73 


6.32 


3 


-0.67 


5.43 


-0.83 


8.35 


4 


-1.10 


14.22 


-1.25 


18.75 


5 


-1.30 


20.36 


-1.46 


25.72 


6 


-1.41 


23.86 


-1.57 


29.63 



Table 3: The coordinates of the individual subarcs 
before (day 0, Figure 1: Top Right) and after the 
shift due to the pulsar proper motion. 
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